/* 
 * File:   VerletMatrix.cpp
 * Author: murmeli
 * 
 * Created on March 18, 2013, 4:15 PM
 */

#include "VerletMatrix.h"

VerletMatrix::VerletMatrix() {
    	setName("Verlet matrix");
        setMatrixMethod(true);
}


VerletMatrix::~VerletMatrix() {
}

void VerletMatrix::methodStep(double& step, vector<double>& xx, vector<double>& yy, vector<double>& dy, cDY* DY, cPot* V) {
    int elements = (int)(V->getUniverse()/step);
    double VStart = V->getUniverseLeft();
    double sq = step*step;
    std::vector<double> tempx;
    //Eigen::LDLT<Eigen::MatrixXd> solver;
    Eigen::MatrixXd mA(elements,elements),mV(elements,elements);
    
    mA.setZero();
    mV.setZero();
 
    tempx.reserve(elements);
    cout << "Filling matrix..." << endl;
 
    mA(0,0) = -2.0/sq;
    mA(0,1) = 1.0/sq;
    
    mV(0,0) = V->getPotentialValueAt(VStart);
    tempx.push_back(V->getPotLeft());
    
    for(int i = 1; i < elements-1; i++ ) {
        mA(i,i-1) = 1.0/sq;
        mA(i,i) = -2.0/sq;
        mA(i,i+1) = 1.0/sq;
        
        mV(i,i) = V->getPotentialValueAt(VStart + i*step);
        tempx.push_back(tempx.back() + step);
    }

    int i = elements - 1;
    mA(i, i-1) = 1.0/sq;
    mA(i,i) = -2.0/sq;
   
    mV(i,i) = V->getPotentialValueAt(VStart + i*step);
    tempx.push_back(tempx.back() + step);
    
    cout << "Done." << endl;
    
    mA = -0.5 * mA + mV;
    setH(mA);
    setX(tempx);
    
}